﻿#include <Geom2d_BSplineCurve.hxx>
#include <Geom_BSplineCurve.hxx>
#include <Geom_BSplineSurface.hxx>
#include <TColStd_Array1OfInteger.hxx>
#include <TColStd_Array1OfReal.hxx>
#include <TColStd_Array2OfReal.hxx>
#include <TColgp_Array1OfPnt.hxx>
#include <TColgp_Array1OfPnt2d.hxx>
#include <TColgp_Array2OfPnt.hxx>
#include <algorithm>
#include <gtest/gtest.h>
#include <iostream>
#include <learn/alg/alg_distance.h>
#include <learn/geom/geom_bspline_curve.h>
#include <learn/geom/geom_bspline_surface.h>
#include <learn/utils/adaptor_occ.h>
#include <learn/utils/output.h>
#include <vector>

using namespace xi::utils;

class TestGeomBSpline : public ::testing::Test
{
  protected:
    std::mt19937 m_random_engine;
    std::uniform_real_distribution<double> m_random_distribution{0, 10};

    void SetUp() override
    {
        m_random_engine.seed(std::random_device{}());
    }

    void TearDown() override{
        // nothing
    };

  public:
    // 分配节点
    void Arrange(const std::vector<double> &knots, TColStd_Array1OfReal &Knots, TColStd_Array1OfInteger &Multiplier)
    {
        Knots.Resize(1, knots.size(), false);
        Multiplier.Resize(1, knots.size(), false);

        Multiplier.Init(0);
        Knots.SetValue(1, knots[0]);
        Multiplier.SetValue(1, 1);

        int j = 1;
        for (int i = 1; i < knots.size(); i++)
        {
            if (fabs(knots[i] - Knots.Value(j)) < 1e-6)
            {
                Multiplier.SetValue(j, Multiplier.Value(j) + 1);
            }
            else
            {
                j++;
                Knots.SetValue(j, knots[i]);
                Multiplier.SetValue(j, 1);
            }
        }
        Knots.Resize(1, j, true);
        Multiplier.Resize(1, j, true);
    }

    // n 个控制点，p 次
    std::tuple<Handle(Geom_BSplineCurve), std::unique_ptr<xi::geom::Geom_BSplineCurve<3>>> GenerateBSplineCurve(int n,
                                                                                                                int p)
    {
        TColgp_Array1OfPnt Poles(1, n);
        std::vector<xi::geom::Vec3> points(n);
        for (int i = 0; i < n; i++)
        {
            points[i] = xi::geom::Vec3(m_random_distribution(m_random_engine), m_random_distribution(m_random_engine),
                                       m_random_distribution(m_random_engine));
            Poles.SetValue(i + 1, gp_Pnt(points[i].x(), points[i].y(), points[i].z()));
        }

        // 构造节点
        std::vector<double> knots;
        for (int i = 0; i < p + 1; i++)
            knots.push_back(0);
        for (int i = p + 1; i < n; i++)
            knots.push_back(knots[i - 1] + m_random_distribution(m_random_engine));

        double final = knots.back() + m_random_distribution(m_random_engine);
        for (int i = n; i < n + p + 1; i++)
            knots.push_back(final);

        TColStd_Array1OfReal Knots;
        TColStd_Array1OfInteger Multiplier;
        Arrange(knots, Knots, Multiplier);

        // 构造曲线
        Handle(Geom_BSplineCurve) occCurve = new Geom_BSplineCurve(Poles, Knots, Multiplier, p);
        std::unique_ptr<xi::geom::Geom_BSplineCurve<3>> curve(
            new xi::geom::Geom_BSplineCurve<3>(std::move(points), std::move(knots), p));
        return {occCurve, std::move(curve)};
    }

    // n 个控制点，p 次
    std::tuple<Handle(Geom_BSplineCurve), std::unique_ptr<xi::geom::Geom_BSplineCurve<3, true>>>
    GenerateBSplineCurve_Weight(int n, int p)
    {
        TColgp_Array1OfPnt Poles(1, n);
        std::vector<xi::geom::Vec4> points(n);
        for (int i = 0; i < n; i++)
        {
            points[i] = xi::geom::Vec4(m_random_distribution(m_random_engine), m_random_distribution(m_random_engine),
                                       m_random_distribution(m_random_engine), m_random_distribution(m_random_engine));
            Poles.SetValue(i + 1, gp_Pnt(points[i].x(), points[i].y(), points[i].z()));
        }

        TColStd_Array1OfReal weights(1, n);
        for (int i = 0; i < n; i++)
        {
            weights.SetValue(i + 1, points[i].w());
        }

        // 构造节点
        std::vector<double> knots;
        for (int i = 0; i < p + 1; i++)
            knots.push_back(0);
        for (int i = p + 1; i < n; i++)
            knots.push_back(knots[i - 1] + m_random_distribution(m_random_engine));

        double final = knots.back() + m_random_distribution(m_random_engine);
        for (int i = n; i < n + p + 1; i++)
            knots.push_back(final);

        TColStd_Array1OfReal Knots;
        TColStd_Array1OfInteger Multiplier;
        Arrange(knots, Knots, Multiplier);

        // 构造曲线
        Handle(Geom_BSplineCurve) occCurve = new Geom_BSplineCurve(Poles, weights, Knots, Multiplier, p);
        std::unique_ptr<xi::geom::Geom_BSplineCurve<3, true>> curve(
            new xi::geom::Geom_BSplineCurve<3, true>(std::move(points), std::move(knots), p));
        return {occCurve, std::move(curve)};
    }

    std::tuple<Handle(Geom_BSplineSurface), std::unique_ptr<xi::geom::Geom_BSplineSurface<3>>> GenerateBSplineSurface(
        int n, int m, int p, int q)
    {
        TColgp_Array2OfPnt Poles(1, n, 1, m);
        std::vector<xi::geom::Vec3> points(n * m);
        for (int i = 0; i < n; i++)
        {
            for (int j = 0; j < m; j++)
            {
                int index = i * m + j;
                points[index] =
                    xi::geom::Vec3(m_random_distribution(m_random_engine), m_random_distribution(m_random_engine),
                                   m_random_distribution(m_random_engine));
                Poles.SetValue(i + 1, j + 1, gp_Pnt(points[index].x(), points[index].y(), points[index].z()));
            }
        }

        // 构造节点
        std::vector<double> uknots;
        for (int i = 0; i < p + 1; i++)
            uknots.push_back(0);
        for (int i = p + 1; i < n; i++)
            uknots.push_back(uknots[i - 1] + m_random_distribution(m_random_engine));

        double ufinal = uknots.back() + m_random_distribution(m_random_engine);
        for (int i = n; i < n + p + 1; i++)
            uknots.push_back(ufinal);

        TColStd_Array1OfReal UKnots;
        TColStd_Array1OfInteger UMultiplier;
        Arrange(uknots, UKnots, UMultiplier);

        std::vector<double> vknots;
        for (int i = 0; i < q + 1; i++)
            vknots.push_back(0);
        for (int i = q + 1; i < m; i++)
            vknots.push_back(vknots[i - 1] + m_random_distribution(m_random_engine));

        double vfinal = vknots.back() + m_random_distribution(m_random_engine);
        for (int i = m; i < m + q + 1; i++)
            vknots.push_back(vfinal);

        TColStd_Array1OfReal VKnots;
        TColStd_Array1OfInteger VMultiplier;
        Arrange(vknots, VKnots, VMultiplier);

        Handle(Geom_BSplineSurface) occSurface =
            new Geom_BSplineSurface(Poles, UKnots, VKnots, UMultiplier, VMultiplier, p, q);
        std::unique_ptr<xi::geom::Geom_BSplineSurface<3>> surface(
            new xi::geom::Geom_BSplineSurface<3>(std::move(points), std::move(uknots), std::move(vknots), p, q));
        return {occSurface, std::move(surface)};
    }

    std::tuple<Handle(Geom_BSplineSurface), std::unique_ptr<xi::geom::Geom_BSplineSurface<3, true>>>
    GenerateBSplineSurface_Weight(int n, int m, int p, int q)
    {
        TColgp_Array2OfPnt Poles(1, n, 1, m);
        std::vector<xi::geom::Vec4> points(n * m);
        for (int i = 0; i < n; i++)
        {
            for (int j = 0; j < m; j++)
            {
                int index = i * m + j;
                points[index] =
                    xi::geom::Vec4(m_random_distribution(m_random_engine), m_random_distribution(m_random_engine),
                                   m_random_distribution(m_random_engine), m_random_distribution(m_random_engine));
                Poles.SetValue(i + 1, j + 1, gp_Pnt(points[index].x(), points[index].y(), points[index].z()));
            }
        }

        TColStd_Array2OfReal weights(1, n, 1, m);
        for (int i = 0; i < n; i++)
        {
            for (int j = 0; j < m; j++)
            {
                int index = i * m + j;
                weights.SetValue(i + 1, j + 1, points[index].w());
            }
        }

        // 构造节点
        std::vector<double> uknots;
        for (int i = 0; i < p + 1; i++)
            uknots.push_back(0);
        for (int i = p + 1; i < n; i++)
            uknots.push_back(uknots[i - 1] + m_random_distribution(m_random_engine));

        double ufinal = uknots.back() + m_random_distribution(m_random_engine);
        for (int i = n; i < n + p + 1; i++)
            uknots.push_back(ufinal);

        TColStd_Array1OfReal UKnots;
        TColStd_Array1OfInteger UMultiplier;
        Arrange(uknots, UKnots, UMultiplier);

        std::vector<double> vknots;
        for (int i = 0; i < q + 1; i++)
            vknots.push_back(0);
        for (int i = q + 1; i < m; i++)
            vknots.push_back(vknots[i - 1] + m_random_distribution(m_random_engine));

        double vfinal = vknots.back() + m_random_distribution(m_random_engine);
        for (int i = m; i < m + q + 1; i++)
            vknots.push_back(vfinal);

        TColStd_Array1OfReal VKnots;
        TColStd_Array1OfInteger VMultiplier;
        Arrange(vknots, VKnots, VMultiplier);

        Handle(Geom_BSplineSurface) occSurface =
            new Geom_BSplineSurface(Poles, weights, UKnots, VKnots, UMultiplier, VMultiplier, p, q);
        std::unique_ptr<xi::geom::Geom_BSplineSurface<3, true>> surface(
            new xi::geom::Geom_BSplineSurface<3, true>(std::move(points), std::move(uknots), std::move(vknots), p, q));
        return {occSurface, std::move(surface)};
    }
};

TEST_F(TestGeomBSpline, curve)
{
    {
        int N = 10;
        int p = 3;
        auto [occCurve, curve] = GenerateBSplineCurve(N, p);

        int M = 1000;
        double first = occCurve->FirstParameter();
        double last = occCurve->LastParameter();
        for (int i = 0; i < M; i++)
        {
            double t = first + (last - first) * i / M;

            gp_Pnt occP;
            gp_Vec occD1, occD2, occD3;
            occCurve->D3(t, occP, occD1, occD2, occD3);

            auto [p, d1, d2, d3] = curve->d3(t);
            xi::geom::Vec3 diffP = p - xi::geom::Vec3(occP.X(), occP.Y(), occP.Z());
            xi::geom::Vec3 diffD1 = d1 - xi::geom::Vec3(occD1.X(), occD1.Y(), occD1.Z());
            xi::geom::Vec3 diffD2 = d2 - xi::geom::Vec3(occD2.X(), occD2.Y(), occD2.Z());
            xi::geom::Vec3 diffD3 = d3 - xi::geom::Vec3(occD3.X(), occD3.Y(), occD3.Z());

            EXPECT_LT((curve->value(t) - p).modulus(), 1e-6);
            EXPECT_LT(diffP.modulus(), 1e-6);
            EXPECT_LT(diffD1.modulus(), 1e-6);
            EXPECT_LT(diffD2.modulus(), 1e-6);
            EXPECT_LT(diffD3.modulus(), 1e-6);
        }

        EXPECT_DOUBLE_EQ(occCurve->FirstParameter(), curve->first_param());
        EXPECT_DOUBLE_EQ(occCurve->LastParameter(), curve->last_param());
        EXPECT_FALSE(curve->Rational);
        EXPECT_EQ(occCurve->Degree(), curve->degree());
        EXPECT_EQ(occCurve->NbPoles(), curve->pole_num());
    }

    {
        int N = 10;
        int p = 3;
        auto [occCurve, curve] = GenerateBSplineCurve_Weight(N, p);

        for (int i = 0; i < curve->pole_num(); i++)
        {
            curve->set_weight(i, m_random_distribution(m_random_engine));
            occCurve->SetWeight(i + 1, curve->weight(i));
        }

        // 升 2 次，降 2 次
        curve->elevate(curve->degree() + 2);
        curve->reduce(1e-6);
        curve->reduce(1e-6);

        int M = 1000;
        double first = occCurve->FirstParameter();
        double last = occCurve->LastParameter();
        for (int i = 0; i < M; i++)
        {
            double t = first + (last - first) * i / M;

            gp_Pnt occP;
            gp_Vec occD1, occD2, occD3;
            occCurve->D3(t, occP, occD1, occD2, occD3);

            auto [p, d1, d2, d3] = curve->d3(t);
            xi::geom::Vec3 diffP = p - xi::geom::Vec3(occP.X(), occP.Y(), occP.Z());
            xi::geom::Vec3 diffD1 = d1 - xi::geom::Vec3(occD1.X(), occD1.Y(), occD1.Z());
            xi::geom::Vec3 diffD2 = d2 - xi::geom::Vec3(occD2.X(), occD2.Y(), occD2.Z());
            xi::geom::Vec3 diffD3 = d3 - xi::geom::Vec3(occD3.X(), occD3.Y(), occD3.Z());

            EXPECT_LT((curve->value(t) - p).modulus(), 1e-6);
            EXPECT_LT(diffP.modulus(), 1e-6);
            EXPECT_LT(diffD1.modulus(), 1e-6);
            EXPECT_LT(diffD2.modulus(), 1e-6);
            EXPECT_LT(diffD3.modulus(), 1e-6);
        }

        EXPECT_DOUBLE_EQ(occCurve->FirstParameter(), curve->first_param());
        EXPECT_DOUBLE_EQ(occCurve->LastParameter(), curve->last_param());
        EXPECT_TRUE(curve->Rational);
    }

    {
        int N = 10;
        int p = 3;
        auto [occCurve, curve] = GenerateBSplineCurve_Weight(N, p);

        int M = 1000;
        double first = occCurve->FirstParameter();
        double last = occCurve->LastParameter();

        // 随机节点插入
        std::uniform_real_distribution<double> dis(first, last);
        std::vector<double> x(10);
        for (int i = 0; i < 10; i++)
        {
            double t = dis(m_random_engine);
            x[i] = t;
            occCurve->InsertKnot(t, 3);
            curve->insert_knot(t, 3);
        }

        // 随机节点细化
        std::vector<double> y(10);
        for (int i = 0; i < 10; i++)
        {
            y[i] = dis(m_random_engine);
            occCurve->InsertKnot(y[i]);
        }
        std::sort(y.begin(), y.end());
        curve->refine_knot(std::move(y));

        // 移除节点
        auto &Mult = occCurve->Multiplicities();
        for (int i = 2; i <= Mult.Length() - 1; i++)
        {
            // occ 的节点删除似乎有点问题
            // occCurve->RemoveKnot(occCurve->Knot(i), Mult[i] - 2, 0);
            curve->remove_knot(occCurve->Knot(i), 2);
        }

        for (int i = 0; i < M; i++)
        {
            double t = first + (last - first) * i / M;

            gp_Pnt occP;
            gp_Vec occD1, occD2, occD3;
            occCurve->D3(t, occP, occD1, occD2, occD3);

            auto [p, d1, d2, d3] = curve->d3(t);
            xi::geom::Vec3 diffP = p - xi::geom::Vec3(occP.X(), occP.Y(), occP.Z());
            xi::geom::Vec3 diffD1 = d1 - xi::geom::Vec3(occD1.X(), occD1.Y(), occD1.Z());
            xi::geom::Vec3 diffD2 = d2 - xi::geom::Vec3(occD2.X(), occD2.Y(), occD2.Z());
            xi::geom::Vec3 diffD3 = d3 - xi::geom::Vec3(occD3.X(), occD3.Y(), occD3.Z());

            EXPECT_LT((curve->value(t) - p).modulus(), 1e-6);
            EXPECT_LT(diffP.modulus(), 1e-6);
            EXPECT_LT(diffD1.modulus(), 1e-6);
            EXPECT_LT(diffD2.modulus(), 1e-6);
            EXPECT_LT(diffD3.modulus(), 1e-6);
        }
    }

    // 测试线性重参数化
    {
        int N = 10;
        int p = 3;
        auto [occCurve, curve] = GenerateBSplineCurve(N, p);

        double first2 = -1;
        double last2 = 1;
        curve->reparam(first2, last2);
        double first = occCurve->FirstParameter();
        double last = occCurve->LastParameter();

        int M = 100;
        for (int i = 0; i < M; i++)
        {
            double t = first + (last - first) * i / M;
            double t2 = first2 + (last2 - first2) * i / M;
            xi::geom::Vec3 p = curve->value(t2);

            gp_Pnt occP;
            occCurve->D0(t, occP);

            auto diffP = p - xi::geom::Vec3(occP.X(), occP.Y(), occP.Z());
            EXPECT_LT(diffP.modulus(), 1e-6);
        }
    }

    // 测试有理重参数化
    {
        std::vector<xi::geom::Vec4> points(6);
        points[0] = xi::geom::Vec4{0.0, 0.0, 0.0, 1.0};
        points[1] = xi::geom::Vec4{1.0, 0.0, 0.0, 0.5};
        points[2] = xi::geom::Vec4{1.0, 1.0, 0.0, 2.0};
        points[3] = xi::geom::Vec4{0.0, 1.0, 0.0, 3.0};
        points[4] = xi::geom::Vec4{0.0, 0.0, 1.0, 0.5};
        points[5] = xi::geom::Vec4{1.0, 0.0, 1.0, 1.0};
        std::vector<double> knots = {0, 0, 0, 0, 0.3, 0.7, 1.0, 1.0, 1.0, 1.0};
        auto points2 = points;
        auto knots2 = knots;
        auto curve = std::make_shared<xi::geom::Geom_BSplineCurve<3, true>>(std::move(points), std::move(knots), 3);
        auto curve2 = std::make_shared<xi::geom::Geom_BSplineCurve<3, true>>(std::move(points2), std::move(knots2), 3);
        curve->reparam(2, 1, 3, 2);

        // std::cout << curve->poles() << std::endl;
        curve->reparam_to_same_weights(1, 2);
        // std::cout << curve->poles() << std::endl;

        std::cout << "distance: " << xi::alg::Alg_Distance<3>::distance(curve, curve2) << std::endl;
    }

    // 测试反转
    {
        int N = 10;
        int p = 3;
        auto [occCurve, curve] = GenerateBSplineCurve(N, p);

        curve->reverse();
        curve->elevate(curve->degree() + 2);
        double first = occCurve->FirstParameter();
        double last = occCurve->LastParameter();

        double first2 = curve->first_param();
        double last2 = curve->last_param();

        int M = 100;
        for (int i = 0; i < M; i++)
        {
            double t = first + (last - first) * i / M;
            double t2 = last2 - (last2 - first2) * i / M;
            xi::geom::Vec3 p = curve->value(t2);

            gp_Pnt occP;
            occCurve->D0(t, occP);

            auto diffP = p - xi::geom::Vec3(occP.X(), occP.Y(), occP.Z());
            EXPECT_LT(diffP.modulus(), 1e-6);
        }
    }
}

TEST_F(TestGeomBSpline, circle)
{
    auto curve = xi::geom::Geom_QR_BSplineCurve<3>::make_circle(
        xi::geom::Vec3{0.0, 0.0, 0.0}, xi::geom::Vec3{0.0, 0.0, 1.0}, xi::geom::Vec3{1.0, 0.0, 0.0}, 10, 45.0, 405.0);
}

TEST_F(TestGeomBSpline, surface)
{
    int Nu = 10, Nv = 8;
    for (int i = 0; i < 4; i++)
    {
        auto [occSurface, surface] = GenerateBSplineSurface(Nu, Nv, 3, 2);

        for (int i = 0; i < Nu; i++)
        {
            surface->pole(i, 4) =
                xi::geom::Vec3(m_random_distribution(m_random_engine), m_random_distribution(m_random_engine),
                               m_random_distribution(m_random_engine));
            occSurface->SetPole(i + 1, 5,
                                gp_Pnt(surface->pole(i, 4).x(), surface->pole(i, 4).y(), surface->pole(i, 4).z()));
        }

        auto [U1, U2, V1, V2] = surface->bounds();

        double occU1, occU2, occV1, occV2;
        occSurface->Bounds(occU1, occU2, occV1, occV2);
        EXPECT_DOUBLE_EQ(occU1, U1);
        EXPECT_DOUBLE_EQ(occU2, U2);
        EXPECT_DOUBLE_EQ(occV1, V1);
        EXPECT_DOUBLE_EQ(occV2, V2);

        int M = 100;
        for (int i = 0; i < M; i++)
        {
            double u = i / (M - 1.0) * (U2 - U1) + U1;
            for (int j = 0; j < M; j++)
            {
                double v = j / (M - 1.0) * (V2 - V1) + V1;

                auto [P, Du, Dv, Duu, Dvv, Duv] = surface->d2(u, v);

                gp_Pnt occP;
                gp_Vec occDu, occDv, occDuu, occDvv, occDuv;
                occSurface->D2(u, v, occP, occDu, occDv, occDuu, occDvv, occDuv);

                xi::geom::Vec3 diffP = P - xi::geom::Vec3(occP.X(), occP.Y(), occP.Z());
                xi::geom::Vec3 diffDu = Du - xi::geom::Vec3(occDu.X(), occDu.Y(), occDu.Z());
                xi::geom::Vec3 diffDv = Dv - xi::geom::Vec3(occDv.X(), occDv.Y(), occDv.Z());
                xi::geom::Vec3 diffDuu = Duu - xi::geom::Vec3(occDuu.X(), occDuu.Y(), occDuu.Z());
                xi::geom::Vec3 diffDvv = Dvv - xi::geom::Vec3(occDvv.X(), occDvv.Y(), occDvv.Z());
                xi::geom::Vec3 diffDuv = Duv - xi::geom::Vec3(occDuv.X(), occDuv.Y(), occDuv.Z());

                EXPECT_LT((surface->value(u, v) - P).modulus(), 1e-6);
                EXPECT_LT(diffP.modulus(), 1e-6);
                EXPECT_LT(diffDu.modulus(), 1e-6);
                EXPECT_LT(diffDv.modulus(), 1e-6);
                EXPECT_LT(diffDuu.modulus(), 1e-6);
                EXPECT_LT(diffDvv.modulus(), 1e-6);
                EXPECT_LT(diffDuv.modulus(), 1e-6);
            }
        }

        EXPECT_FALSE(surface->Rational);
        EXPECT_EQ(occSurface->NbUPoles(), surface->u_pole_num());
        EXPECT_EQ(occSurface->NbVPoles(), surface->v_pole_num());
        EXPECT_EQ(occSurface->UDegree(), surface->u_degree());
        EXPECT_EQ(occSurface->VDegree(), surface->v_degree());

        std::swap(Nu, Nv);
    }

    for (int i = 0; i < 4; i++)
    {
        auto [occSurface, surface] = GenerateBSplineSurface_Weight(Nu, Nv, 3, 2);

        for (int i = 0; i < Nu; i++)
        {
            surface->set_weight(i, 4, m_random_distribution(m_random_engine));
            occSurface->SetWeight(i + 1, 5, surface->weight(i, 4));
        }

        auto [U1, U2, V1, V2] = surface->bounds();

        // 随机节点插入
        std::uniform_real_distribution<double> disU(U1, U2);
        std::uniform_real_distribution<double> disV(V1, V2);
        for (int i = 0; i < 10; i++)
        {
            double u = disU(m_random_engine);
            double v = disV(m_random_engine);
            occSurface->InsertUKnot(u, 2, 0);
            surface->insert_u_knot(u, 2);

            occSurface->InsertVKnot(v, 2, 0);
            surface->insert_v_knot(v, 2);
        }

        // 随机节点细化
        std::vector<double> x(10), y(10);
        for (int i = 0; i < 10; i++)
        {
            x[i] = disU(m_random_engine);
            y[i] = disV(m_random_engine);
            occSurface->InsertUKnot(x[i], 1, 0);
            occSurface->InsertVKnot(y[i], 1, 0);
        }
        std::sort(x.begin(), x.end());
        std::sort(y.begin(), y.end());
        surface->refine_u_knot(std::move(x));
        surface->refine_v_knot(std::move(y));

        double occU1, occU2, occV1, occV2;
        occSurface->Bounds(occU1, occU2, occV1, occV2);
        EXPECT_DOUBLE_EQ(occU1, U1);
        EXPECT_DOUBLE_EQ(occU2, U2);
        EXPECT_DOUBLE_EQ(occV1, V1);
        EXPECT_DOUBLE_EQ(occV2, V2);

        int M = 100;
        for (int i = 0; i < M; i++)
        {
            double u = i / (M - 1.0) * (U2 - U1) + U1;
            for (int j = 0; j < M; j++)
            {
                double v = j / (M - 1.0) * (V2 - V1) + V1;

                auto [P, Du, Dv, Duu, Dvv, Duv] = surface->d2(u, v);

                gp_Pnt occP;
                gp_Vec occDu, occDv, occDuu, occDvv, occDuv;
                occSurface->D2(u, v, occP, occDu, occDv, occDuu, occDvv, occDuv);

                xi::geom::Vec3 diffP = P - xi::geom::Vec3(occP.X(), occP.Y(), occP.Z());
                xi::geom::Vec3 diffDu = Du - xi::geom::Vec3(occDu.X(), occDu.Y(), occDu.Z());
                xi::geom::Vec3 diffDv = Dv - xi::geom::Vec3(occDv.X(), occDv.Y(), occDv.Z());
                xi::geom::Vec3 diffDuu = Duu - xi::geom::Vec3(occDuu.X(), occDuu.Y(), occDuu.Z());
                xi::geom::Vec3 diffDvv = Dvv - xi::geom::Vec3(occDvv.X(), occDvv.Y(), occDvv.Z());
                xi::geom::Vec3 diffDuv = Duv - xi::geom::Vec3(occDuv.X(), occDuv.Y(), occDuv.Z());

                EXPECT_LT((surface->value(u, v) - P).modulus(), 1e-6);
                EXPECT_LT(diffP.modulus(), 1e-6);
                EXPECT_LT(diffDu.modulus(), 1e-6);
                EXPECT_LT(diffDv.modulus(), 1e-6);
                EXPECT_LT(diffDuu.modulus(), 1e-6);
                EXPECT_LT(diffDvv.modulus(), 1e-6);
                EXPECT_LT(diffDuv.modulus(), 1e-6);
            }
        }

        EXPECT_TRUE(surface->Rational);
        EXPECT_EQ(occSurface->NbUPoles(), surface->u_pole_num());
        EXPECT_EQ(occSurface->NbVPoles(), surface->v_pole_num());
        EXPECT_EQ(occSurface->UDegree(), surface->u_degree());
        EXPECT_EQ(occSurface->VDegree(), surface->v_degree());

        std::swap(Nu, Nv);
    }

    // 升阶
    for (int i = 0; i < 4; i++)
    {
        auto [occSurface, surface] = GenerateBSplineSurface_Weight(Nu, Nv, 3, 2);

        // 曲面升阶
        occSurface->IncreaseDegree(5, 4);
        surface->elevate(5, 4);

        auto [U1, U2, V1, V2] = surface->bounds();
        double occU1, occU2, occV1, occV2;
        occSurface->Bounds(occU1, occU2, occV1, occV2);
        EXPECT_DOUBLE_EQ(occU1, U1);
        EXPECT_DOUBLE_EQ(occU2, U2);
        EXPECT_DOUBLE_EQ(occV1, V1);
        EXPECT_DOUBLE_EQ(occV2, V2);

        int M = 100;
        for (int i = 0; i < M; i++)
        {
            double u = i / (M - 1.0) * (U2 - U1) + U1;
            for (int j = 0; j < M; j++)
            {
                double v = j / (M - 1.0) * (V2 - V1) + V1;

                auto [P, Du, Dv, Duu, Dvv, Duv] = surface->d2(u, v);

                gp_Pnt occP;
                gp_Vec occDu, occDv, occDuu, occDvv, occDuv;
                occSurface->D2(u, v, occP, occDu, occDv, occDuu, occDvv, occDuv);

                xi::geom::Vec3 diffP = P - xi::geom::Vec3(occP.X(), occP.Y(), occP.Z());
                xi::geom::Vec3 diffDu = Du - xi::geom::Vec3(occDu.X(), occDu.Y(), occDu.Z());
                xi::geom::Vec3 diffDv = Dv - xi::geom::Vec3(occDv.X(), occDv.Y(), occDv.Z());
                xi::geom::Vec3 diffDuu = Duu - xi::geom::Vec3(occDuu.X(), occDuu.Y(), occDuu.Z());
                xi::geom::Vec3 diffDvv = Dvv - xi::geom::Vec3(occDvv.X(), occDvv.Y(), occDvv.Z());
                xi::geom::Vec3 diffDuv = Duv - xi::geom::Vec3(occDuv.X(), occDuv.Y(), occDuv.Z());

                EXPECT_LT((surface->value(u, v) - P).modulus(), 1e-6);
                EXPECT_LT(diffP.modulus(), 1e-6);
                EXPECT_LT(diffDu.modulus(), 1e-6);
                EXPECT_LT(diffDv.modulus(), 1e-6);
                EXPECT_LT(diffDuu.modulus(), 1e-6);
                EXPECT_LT(diffDvv.modulus(), 1e-6);
                EXPECT_LT(diffDuv.modulus(), 1e-6);
            }
        }

        EXPECT_TRUE(surface->Rational);
        EXPECT_EQ(occSurface->NbUPoles(), surface->u_pole_num());
        EXPECT_EQ(occSurface->NbVPoles(), surface->v_pole_num());
        EXPECT_EQ(occSurface->UDegree(), surface->u_degree());
        EXPECT_EQ(occSurface->VDegree(), surface->v_degree());

        std::swap(Nu, Nv);
    }
}